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SECTION  I 


INTRODUCTION 

For  some  lightly  damped  structure  the  free  response  may  be  expressed  as  a  sum  of  a 
number  of  complex  exponentials  with  complex  amplitudes.  It  is  clear  that  the  number  of 
such  terms  is  not  necessarily  the  same  for  different  experiments  on  the  same  structure. 
The  number  of  these  terms  depends  upon  the  initial  conditions.  From  a  practical  point  ol 
view  one  wxmld  be  satisfied,  usually,  with  a  free  response  which  contains  all  the  frequencies 
in  some  limited  range  of  interest. 

In  Section  II  we  give  the  details  for  efficiently  and  accurately  carrying  out  Prony’s 
method,  [l],  for  determining  the  complex  exponentials  and  complex  amplitudes  from  evenly 
spaced  sampled  data  assuming  that  the  number  of  such  terms  present  is  known.  The 
method  given  is  also  applicable  when  the  data  is  averaged  in  an  appropriate  but  natural 
way.  There  is  no  difficulty  in  applying  the  method  to  multidegree  of  freedom  cases.  We 
believe  that  one  could  possibly  treat  the  case  of  a  complex  frequency  wuth  multiplicity 
greater  than  1. 

In  Section  III  a  procedure,  the  Progressive  Algorithm  [2],  is  given  which  essentially 
enables  us  to  determine  the  number  of  complex  frequencies  present.  Actually  the  progres¬ 
sive  algorithm  performs  the  first  step  of  the  method  given  in  Section  II.  Thus,  it  is  no 
longer  necessary  to  know  in  advance  the  number  of  complex  frequencies  present.  The  pro¬ 
gressive  algorithm  is  also  applicable  to  the  averaged  data.  Given  a  real  valued  function  of 
a  finite  number  of  distinct  complex  frequencies  w  ith  their  associated  complex  amplitudes 
it  is  clear  that  a  method  comprised  of  the  procedures  given  in  Sections  II  and  III  is  a 
theoretically  exact  method  for  determining  these  frequencies  and  amplitudes  from  equally 
spaced  function  values. 

In  Section  IV  w'e  derive  the  algebraic  properties  associated  with  the  method  of  Least 
Squares  for  the  case  of  one  quantity  depending  linearly  on  another.  This  in  turn  provides 


the  background  for  developing  the  method  of  Lattice  Filters.  [3).  for  stationary  sequences 
and  the  method  of  Canonical  Variate  Analysis.  [4]. 
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SECTION  II 

DETERMINATION  OF  COMPLEX  EXPONENTIALS 

Consider  a  sum  of  terms  of  the  form  a*e;rp(A*/).  The  coefficient  may  be  either  a 
vector  or  scalar.  For  every  such  term  there  is  a  companion  term  which  is  the  complex  con¬ 
jugate.  so  that  the  components  of  the  sum  are  real.  Such  a  sum  arises  as  the  homogeneous 
solution  of  a  linear  system  of  ordinary  differential  equations  with  constant  coefficients.  In 
particular  dynamical  systems  with  symmetric  mass,  stiffnes  and  damping  matrices.  Re¬ 
gardless  of  origin,  the  problem  of  determining  the  parameters  a*  and  A*  from  experimental 
data  is  frequently  encountered. 

In  this  section  we  examine  and  discuss  a  very  old  method  associated  with  the  name 
Prony.  The  method  or  procedure  described  is  efficient,  easily  programmed  and  very  ac¬ 
curate  when  tested  on  numerically  produced  data  sampled  at  an  appropriate  rate.  In  a 
practical  problem  the  number  of  the  terms  a*e.rp(A*/)  is  not  known.  In  the  next  section 
we  give  a  method  for  dealing  with  this  difficulty.  There  is  the  problem  of  multiplicity  for 
a  particular  A*  which  we  will  discuss  also. 

Consider  a  real  valued  function  x(t)  defined  by 

fl 

■r(0  =  Y  a*exp(A*0  (1) 

k  1 

where  n  =  2m, 

am+k  —  Q k 

and 

=  Xf;. 

The  bar  denotes  the  complex  conjugate.  The  Eq  (l)  can  be  written  also  as 
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»/Ar  =  exp(Xkh) 


From  Eq  (1)  wo  obtain 


s(i)  =  die  '  +  •••  +  a„/*' 
x(l  +  h)  =  r/,o1eA''  +  -  ••  +  ilna„rXr't 

x(t  +  nh)  —  + - h  r)ZaneXnt. 

The  Eqs  (4)  may  be  written  as 


x(t) 

-1 

.  -i 

x(t  +  h) 

-m 

•  -»?n 

x(t  4-  nh) 

-*r 

•  -nS 

Since  this  homogeneous  system  of  equations  has  a  nontrivial  solution  the  determinant  of 
the  system  matrix  is  zero. 

Let  Ajk  denote  the  cofactor  of  the  jkLh  element  in  the  determinant  of  the  system 
matrix.  Eq  (5).  \Ne  have 

-4„^i  i x(t  +  nh)  +  A„  , x(t  +  (n  -  1  )h)  +  •  ••  +  A,  ,x(f)  -  0.  (6) 

Assume  for  the  present  that  the  A*  and  hence  also  the  r)k  are  different  from  one  another. 
Observe.  Eq  (5).  that  the  minor  of  the  element  x(t  +  nh)  is  a  Vandermonde  determinant. 
Hence  the  cofactor  <4„_  1 1  ^  0  and  Eq  (6)  can  be  written  as 

x{t  +  nh)  +  +  (n  -  1)/?)  + - (-  -11  j(/)  =  0.  (7) 

The  Eqs  (6)  and  (7)  show  that  the  function  values  x(l).  x(t  +  h).  ...  satisfy  a  linear  nth 
order  difference  relation. 

In  Eq  (7)  set  the  coefficient 
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and  rewrite  Eq  (7)  as 


x(l  +  nh)  +  p„x(t  +  (n  -  1  )//)+•••  +  ;i,x(/)  =  0. 


In  Eq  (5)  replace  the  first  column  of  the  coefficient  matrix  by  the  symbols 


1 .  —  —  t/2.  ■  ■  ■ .  —Tjn 


Then  in  the  same  way  that  we  obtained  Eq  (9).  the  expansion  of  the  determinant  of  the 
modified  coefficient  matrix  will  give  the  polynomial 


V”  +  PnVn  1  +  •  •  ■  +  P\ 


and  it  is  clear  from  elementary  determinant  theory  that  the  r/*,  Eq  (3).  are  the  roots  of 
the  polynomial  equation 

P{r>)  =  P"  +  Pn1n -'  +  *■•  +  Pi  =0  (10) 

Nou  from  elementary  college  algebra  we  have 

P{v)  =  n(l- 

k-l 

m 

=  ~  ***) 
k  1 

=  n  (>?2  -  2ffPeM + (PeM2 + (imirik})2) 

k-  1 

since  the  n  roots  of  Eq  (10)  occur  in  complex  conjugate  pairs.  From  this  it  is  clear  that 
the  coefficients.  /»* .  of  Eq  (9).  and  also  in  Eq  (10).  are  real. 

Set  x(i)  =  X],  r(t  +  h)  =  x2 . x(t  +  nh)  =  xn^\.  and  so  on.  From  Eq  (9)  we  obtain 


the  system  of  equations 
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Thus,  in  principle,  from  a  set  of  2n  values  of  the  function  x(t)  the  coefficients  />,.  ;»2 . 

l>„  ran  be  determined,  if  the  coefficient  matrix  in  Eq  (11)  is  nonsingular.  From  Eq  (4)  we 
obtain  ,  , 


and  similarly 
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That  is. 

we  obtain 
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Now  the  determinant  of  the  left  hand  side  is  equal  to  the  product  of  the  determinants  of 

the  right  hand  side.  By  our  assumptions  above,  it  is  clear  that  the  value  of  each  of  the 
determinants  on  the  right  hand  side  is  different  from  zero.  It  follows  from  our  observations. 
Eq  (12),  that  we  may  assume  that  the  coefficients  are  determined  by  the  Eqs  (11)  or 
known. 

We  consider  next  the  second  task,  namely  determining  the  zeros  of  the  polynomial 
P(>])-  The  zeros  r/*  of  the  polynomial  P{rj)  occur  in  complex  conjugate  pairs.  Hence  we 
only  need  to  determine  the  zeros  lying  in  the  upper  half  plane  For  A*.  =  o*  +  i 3k  we  have 
from  Eq  (3) 

r,k  =  r,*»V,*h.  (13) 
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In  the  present  discussion  we  are  only  considering  lightly  damped  systems.  Thus  ak  <  0 
and  is  numerically  small.  Hence  0  <  rxp(nkh)  <  1  and  it  is  clear  from  Eq  (13)  that  all  the 
zeros  T]k  lie  inside  the  unit  circle.  In  particular,  if  nk  is  sufficiently  small  numerically  and 
and  if  h  is  sufficiently  small  also  then  the  zeros  i)k  will  he  close  to  1  in  magnitude.  That 
is  the  roots  r/k  lie  close  to  the  circumference  of  the  unit  circle. 

If  the  zeros  rjk  lie  close  to  the  unit  circle  we  are  able  to  determine  them  to  a  high  degree 
of  accuracy  in  two  “steps".  The  first  step  is  a  linear  search  along  the  circumference  of  the 
unit  circle  which  gives  an  initial  approximation  f)k  to  the  zero  r)k.  The  second  step  is  a 
Newton  iteration  process  starting  with  this  initial  approximation. 

To  implement  the  first  step  we  have 

1  dP  1  1 

P(v)  dtl  r)  ~  Vk  V  -  rln  ’ 


\dP/drj\ 


!  Pin)  I' 

By  means  of  a  search,  determine  approximately  local  maxima  rjk  of  /(?/)  along  the  cir¬ 
cumference  of  the  unit  circle  in  the  upper  half-plane.  One  should  find  m  values  of  tjk. 
1  <  k  <  m.  It  is  clear  that  the  search  is  not  restricted  to  the  circumference  of  the  unit 
circle.  One  may  search  along  the  circumference  of  circles  of  radii  different  from  1. 

The  recurrence  relation  for  a  Newton  Process  for  determining  an  r/  satisfying  /’(q)  =  0 


^  1  r> 


P(h) 
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z0  -  *lk 


where  P'(r))  denotes  the  derivative  of  P(t j)  with  respect  to  q. 
Let  rfh  denote  one  of  these  roots.  We  have 
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„  _  _Ai/j 

rfk  =  € 


which  we  write  as 


=  r _  eX„h  _  e{ok  -.,<*)/> 


Vk  =  rk 


From  this  we  obtain 


i|.  it.  41. 


(H  =  (l/h)/oy(r*  ] 


^  =  («*  +  2; r/)//,. 


It  is  clear  from  Eq  (18)  that  A*,  is  not  uniquely  defined.  One  takes  /  =  0  in  Eq  (18)  unless 


tliere  is  good  reason  to  do  otherwise. 


Once  the  A*.'s  are  determined,  it.  only  remains  to  determine  the  o*'s.  The  followii 


system  of  equations  is  available  for  determining  the  ak  s 


V\  V2 


1  a1  e 

T)n  a2ex 


This  system  can  be  solved  by  a  linear  systems  solver.  However,  because  of  the  special 
form  of  the  coefficient  matrix  one  can  determine  the  t'alue  of  the  product  akeXkt  in  a  more 


efficient  wav. 


Here  now  we  develop  the  procedure  for  solving  the  system,  Eqs  (19),  which  utilizes 
the  special  form  of  the  coefficient  matrix.  Denote  the  cofactor  of  the  j  kth  element  of  the 


coefficient  matrix  by  B}k.  Then  by  Cramer's  rule 


(Yi  <  ' B)k)okex =  £  x}D}k 
j-t  j:  i 


}~  1 

Divide  Eq  (20)  through  by  Bnk.  The  coefficient  of  akex,',  in  Eq  (20)  is  a  polynomial  in  r/ 


of  degree  n  -  1  evaluated  at  t]  =  r] k.  We  have 


pi  „  \  _  „n  1  .  -  1 k  n  2  ,  ,  k 


-  y  8j*  j  i 

^  B„  k  ?/*  ■ 


It  is  clear  from  Eq  (19)  that  the  roots  of  the  polynomial  equation  P(i))  -  0  are  the  r]} , 


j  ^  k.  Hence 


r(’i)  =  n  -  rb )• 


The  polynomial,  Eq  (10),  is 


P{v)  —  v"  +  Pn1/"  +  •  ■  •  +  ;>) 

n 

=  11^  -  »o>- 

j-  i 


Then 


—  =  nr?"  +  (r?  -  1  )/>„!?"  +  •  -  •  +  p2 

t=i  j=i 

=  nr?”-1  +  (n  -  l)pnv”~2  +  •••  +  p2 

df]  '»J  =  «7*  K  1 

=  n  (»?*  -  *?/)  =  PM- 


Thus  the  coefficient  of  ate-***  is 


x  <*P; 

p(,‘>  = 


The  right  hand  side,  RHS,  of  Eq  (20)  after  division  by  is  the  quantity 


RHS=  Yx(t  +  (j  -  1)/>)|J- 

3  1 


That  is,  /?#5  is  the  polynomial  />(?7)  with  r\3  replaced  by  the  value  x(t  +  jh). where 
j  =  0, - n  —  1.  We  make  the  following  observations: 

P(»?)  =  P(n)/{r)  ~  *1k) 

=  s5nce  PM  =  0  (28) 

=  ^^+Pn^.i_!?r!  +  ...+P2. 

»?-»?*  »?  -  Vk 

For  the  case  n  =  4  Eq  (28)  becomes 

P(v)  =  r)3  +  V2rik  +  *?*?*  +  Vk 
+  Pa{v2  +  Wk  +  vik) 

(29) 

+  P3(h  +  hfc) 
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Now  we  rewrite  P(t/ )  given  by  Eq  (29)  as  a  polynomial  in  r/*.  We  have 

/>(»>)  -  >n 

■+  (n  +  Pa  )  w 

(3°) 

+  (>l2  +  PaP  +  Ps)Vk 
+  (v3  +  PaV2  +  P:\P  +  Pi)- 

Finally  to  obtain  the  value  which  we  called  RHS  replace  r/J  in  Eq  (30)  by  r(/  +  jh),  where 

j  —  0 . n  —  1  and  obtain,  for  the  case  n  =  4, 

RHS  =  r(t)4 

+  (rii  +  h)  +  p4x(t))r)l 

(31) 

+  (*(/  +  2 h)  +  p4x(t  +  h)  +  P3x[t))  Pic 
+  x(t  +  3/j)  +  p4x(t  +  2h)  -f  p3i(<  +  h)  +  p2x(t). 

The  expression  for  the  value  of  RHS  for  a  general  value  of  n  is  readily  obtained  from  Eq 


(31).  We  have 


RHS  =  J-(/)^n 


+  (j(/  +  //)  +  p„x(t))i}"  2 
+  (r(t  +  2 h)  +  p„x(l  +  fe)  +  pn.  ,x(f))ttf 

-  x(t  +  (n  -  l)/i)  +  p„x(t  +  (n  -  2 )h)  +  •••-+  p2r{t). 


Now  we  have  shown  that 


is  equal  to  the  quotient  of  two  polynomials  of  degree  n  —  1.  each  of  which  is  evaluated  at 
r;*. .  Polynomials  are  readily  and  rapidly  evaluated.  The  numerator  is  given  by  Eq  (32)  and 
the  denominator  by  Eq  (26).  The  polynomials  are  the  same,  have  the  same  coefficients, 
for  all  values  of  k.  for  1  <  k  <  m  and.  of  course. 


nk  .  m  -  °k- 


I'sually  there  is  no  reason  for  not  taking  1  —  0  in  Eq  (33). 

There  is  no  difficulty  in  going  from  the  scalar  case  just  discussed  to  the  case  where 
the  coefficients  a*  are  p-dimensional  vectors.  One  applies  the  procedure  described  above 
to  a  component.  The  A*  do  not  have  to  be  determined  for  each  component.  Thus  the 
determination  of  the  polynomial  P{tj)  is  done  for  only  one  component.  The  coefficients  for 
the  polynomial.  Eq  (32).  have  to  be  determined  for  each  component  using  the  observed 
values  for  that  component. 

Reviewing  the  above  development  of  a  procedure  for  the  determination  of  complex 
exponentials  we  see  that  the  first  task  is  the  determination  of  the  coefficients  pk  of  a 
difference  equation,  Eq  (9).  These  coefficients  are  also  the  coefficients  for  a  polynomial 
equation,  Eq  (10).  The  second  task  is  the  determination  of  the  roots  r)k  of  this  polynomial 
equation.  Once  the  r)k  are  known  the  A*  are  also  determined  by  Eqs  (17)  and  (18).  The 
last  task  is  the  determination  of  the  ak's  or  the  a*'s  using  Eqs  (26)  and  (32).  If  the  number 
n  of  complex  exponentials  is  known  the  first  task  can  be  accomplished  by  solving  a  system 
of  linear  equations.  Eq  (11).  For  the  case  n  not  known  an  alternative  to  the  use  of  Eq  (11) 
is  given  in  SECTION  III. 

For  the  vector  case  there  is  the  possibility  of  two  or  more  linearly  independent  vectors 
a*  as  the  coefficient  of  the  same  complex  exponential  in  the  true  expression  for  the  function 
f(/).  Eq  (1).  It  is  clear,  however,  that  the  experimental  data  can  only  express  the  sum. 
that  is,  these  two  or  more  coefficients  of  the  same  complex  exponential  are  summed  and 
regarded  as  one.  If  one  suspects  that  there  are  multiple  roots  one  must  process  other  sets  of 
data  with  different,  initial  conditions.  Let  a*j  and  a*2  denote  the  coefficients  corresponding 
to  the  value  A*  obtained  from  two  different  sets  of  data.  If  A*  is  not  a  multiple  root  a*2 
will  be  some  scalar  multiple  of  a*j.  On  the  other  hand  if  A*  is  a  multiple  root  then 

«fc2  ~  <*«l H  /  9 

for  any  value  of  the  scalar  c.  Hence  one  must  process  a  number  of  data  sets  equal  to  the 
multiplicity  of  the  A*. 
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We  have  shown  that  a  set  of  2v  equally  spared  values  of  the  function  j-(/),  Eq  (1), 
completely  determine  the  function  x(t)  in  principle.  Next  we  show  that  the  algorithms 
developed  for  determining  the  function  x(/)  apply  also  to  appropriate  average $  of  the 
evenly  spared  values  of  x(/). 

Suppose  that  we  have  very  many,  say  JV  >>  2 n.  equally  spaced  values  of  the  function 
x[t).  From  Eq  (11)  we  have,  for  example 

x,p,  +  x2p2  +  •  •  +  xnPn  =  -t,h  , 

X2P\  +  x'iP2  +  ‘  ‘  '  +  Jn-i|  —  ~xn\2 

xtP\  "F  xr  *  lP2  *F  '  ’  ‘  "F  ♦  r  tPn  =  ~xn  4  r- 
Adding  these  equations  together  and  dividing  by  r  we  obtain 

1  r,  1  r,  1  !  r 

(r  Y  *k)P\  +  (r  Y  xk*\)P2  +  •••  +  (r  Y  1  )Pn  ^  -  rYTn-k 

k  1  k  I  k~ 1  k  I 

S('t 

1  r 

*i  =  r  Y  r*  * 3  1-  (34) 

k  1 

Thus  replacing  the  value  x}  in  Eq  (11)  by  the  average  x}  we  obtain  equations  of  exactly 
the  same  form  as  Eqs  (11)  for  determining  the  coefficients  pk. 

From  Eq  (4).  see  also  the  equations  immediately  preceding  and  leading  up  to  Eq  (12). 
we  have 


x\ 

xi 

xr 

x2 

x  > 

Xr>, 

xn 

xn •*  1  •  ■ 

xn~r j 
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Pi 

P2 

Pn 

n” 

'  P2"  ‘ 

•••  Pr7 

0|CAl' 

njfA'fp,  . 

a]ex',tf] 

xi  x 2  ^ 

a2rx?',h  . 

..  «2eA,,p; 

ar)rA”f 

<‘„<K,Pn 

aneXy'i  i)T„ 

Multiplying  this  equation,  both  sides,  on  the  right  by  an  r  vector  with  all  components 
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equal  to  1  and  dividing  through  by  r  we  obtain 


m  m 


Li-„J  L  T]"~  1  r;J  1  ...  >>"  'j  L  r  (1  +  »?«  +  •••  +  >>rV  V»V"'J 

Thus  when  using  average  values  as  given  by  Eq  (34)  we  obtain,  instead  of  Eq  (19)  which 
wo  solved  for  the  product  ak^k> .  Eq  (36)  which  is  solved  for  the  quantity 

*(1  +  m  +  •■•  +  ’ll  1  )oicCXl,t  ■ 

The  algorithm  given  for  solving  for  the  value  of  is  also  applicable  for  solving  for  the 
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SECTION  III 

THE  PROGRESSIVE  ALGORITHM 


In  this  section  we  consider  the  progressive  algorithm  which  was  used  by  C.  Lanczos  in 
[2j.  This  procedure  is  useful  when  the  number  of  complex  exponentials  in  the  function,  Eq 
(1).  is  not  known.  To  develop  the  algorithm  we  first  describe  the  procedure  from  the  point 
of  view  of  systems  of  equations.  Then  we  repeat  the  description  using  matrix  notation. 
The  essential  computations  are  summarized  in  Eqs  (55-60). 

For  discussion  purposes  we  assume  a  set  of  numbers  x *  and  suppose  there  is  a  least 
value  n  and  a  vector  having  components  »?i, . . . ,  1  satisfying  the  condition 


x\ 

*2 

1 2 

X 3 

xn 

xn  4  1 

xn  i  1 

■r2n  I 


Such  a  system  of  homogeneous  equations  arises  in  the  obvious  way  from  Eq  (9).  Typically 
for  a  problem  of  this  nature  one  might  delete  a  row.  say  the  last  row'  of  this  system,  Eq 
(37).  take  the  entries  remaining  in  the  last  column  to  the  other  side  of  the  equal  sign  and 
in  this  way  obtain  the  system  of  equations,  Eq  (11).  The  difficulty  is  that  the  value  of  n 
is  not  knowm.  We  could  still  use  such  an  attack  by  taking,  in  turn,  n  =  3.4. .. .  until  we 
obtained  a  value  of  n  for  which  Eq  (37)  was  satisfied  in  some  numerical  sense.  The  special 
form  of  the  coefficient  matrix.  Eq  (37),  enables  one  to  develop  an  efficient  procedure  for 
carrying  out  this  approach. 

Suppose  wre  have  obtained,  some  how.  for  a  value  of  k  <  n  the  following  conditions 


J’l£l+J'2£2  +  ',‘  +  -r*£/t  +  -r*:4l  -  0 

•r2^l  +-r3s2  +  -  -  +  JrJt-ls<r+Jt-2-0 


+  xk  -  2*2  +  '  ’  ‘  ■*-  x2k£k  +  -r2/t-l  —  ^1 


* 

i . 
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x2<i\  +  x3?2  +  •  ‘  +  xk  +  2  —  0 

■T3?l  +  x4$2  +  '  •  •  +  xk  +  2$k  +  xk-t3  ~  0 


xk  +  2?1  +  xk  +  3$2  +  •  *  ’  +  x2k^  1  Clr  +  x2k^2  —  ^2 


Observe  that  by  multiplying  each  of  the  Eqs  (38)  by  -h2/hx  and  adding  the  resulting 
system  of  equations  to  the  Eqs  (39)  one  obtains 


X]V\  +  x2rl2  +  •••  +  Tk- \Vk  >,  +  xk+ 2  =  0 
x2rh  +  x3rl2  +  '  '  •  +  x k 2r)k~i  \  +  xk  +  3  =  0 


•ri  +  l,7l  +  xk-+2rh  +  ‘  ‘  '  +  x2k+  l7?*-*  1  +  I2*  +  2  -  0 


where 


Next  compute 


V\  _  (~h2/h j)£i 
V2  —  (~^2/h\)i2  +  ft 


Vk  -  {-h2/hi)(,k  +  f*-  i 
*?*+!  -  (  ^2/ ^  1 )  +  ft- 


xk  +  2rl\  +  xk4  3rl2  +  ‘  ‘  ‘  +  x2k*  2rlk*  1  +  J2t  +  3  -  h. 


If  h  =  0  then  the  set  of  equations  comprised  of  the  Eqs  (40)  and  Eq  (42)  is  the  system 
Eq  (37)  and  the  desired  solution  vector  is  given  by  the  Eqs  (41).  If  h  f-0  then  the  system. 
Eqs  (40)  and  Eq  (42),  will  play  the  role  of  the  system.  Eq  (38).  We  need  to  obtain  the 
system  which  will  play  the  role  of  the  system.  Eq  (39).  To  obtain  the  new  system  Eq  (39) 
multiply  the  old  Eq  (39)  by  —h/h 2  and  add  to  the  new  system  Eq  (38)  with  the  first  row 
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deleted.  We  have 

x2(V\  -  fiVM  +  *3(*?2  -  ^h/h2)  +  ■■■  +  It,2(»?jt+1  -  hjh2)  +  xic-r 3  =  o 
X3(Vl  -  f)  +  X4(rh  -  ^h/h2)  +  •  •  ■  +  X*.+3(^  +  ,  -  /)//l2)  +  Xt  +  4  =  0 


3*  +  2(»?l  -  ?1^/M  +  **43(»?2  -  ftVM  +  ••• +l2fc  +  2('7Jt-l  -  VM  +  3-2*4  3  =  0. 

From  the  Eqs  (43)  set 

fli  =  Vi  ~  Cifc/^2 
fh  =  m  -  &h/h2 

:  (44) 

nk-rtk-  Skh/h 2 
>?*-n  =  »?*4i  -  V^- 

l'sing  the  values  obtained  in  the  Eqs  (44)  compute 

3*  +  3*)i  +  Zk^f)2  +  •••  +  x2*+3')*4i  +  32*m  =  h.  (45) 

Replacing  the  q's  by  £'.s,  /i  by  ft],  the  i)'s  by  $'s  and  h  by  h2  the  Eqs  (40)  and  (42)  are  of 
the  form  Eqs  (38)  and  Eqs  (43)  and  (45)  are  of  the  form  Eqs  (39).  We  are  now  ready  to 
repeat  the  process  described  above. 

The  Eqs  (38-45)  completely  describe  the  Progressive  Algorithm,  but  from  these  equa¬ 
tions  it  is  somewhat  difficult  to  see  what  is  going  on.  Here  now  we  will  develop  the 
progressive  algorithm  using  matrix  notation.  The  Eqs  (38)  and  (39)  can  be  written  as 


where 


A't,  = 


A*  2  — 


An 

f7*  =  #*, 

A'*  2 

n  =  tf*2 

3-1 

3*2 

3**  r  1 

3-2 

3-3 

3**-  2 

.3*- J 

J-/  -2 

J2*-»  1 

3-  2 

3*3 

3*-2 

3-3 

3  4 

3*  +  3 

.J"* -2 

3**  -  3 

32*-2 

1 

16 

and 


uk  =  ■  ,  Hkl  = 

&  o 

-  1  -  . 


:  .  r*  =  ■  .  Hk  2  = 

o  ft 


We  think  of  the  matrix  A”*!  as  partitioned  into  two  block  submatrices 


■  *1  ■ 

*2 

xk^i  ' 

AT  11  = 

x2 

,  at  12  = 

*3 

xk^2 

■  xk+  1- 

,xk  ~2  ■ 

x2k+\  . 

The  block  submatrices  of  the  matrix  A”*  2  are 


AT  21  = 


x2  ■  ■  ■  xk- f  1 
x3  2-*  +  o 

.  xk  +  2  ■  ■  ■  x2k+  1 


>  Xk  22  — 


.  x2k  +  2 


Observe  that 


Xk  12  -  Xk  21  • 

The  vectors  Uk  and  Vk  are  partitioned  in  the  obvious  way  so  as  to  be  compatible  with  the 
partitioning  of  the  matrices  AT  1  and  AT  2-  We  have 

£2' 

r*i  =  Ui].r*2=  j 

{ft 
1  . 

and 

ft 

^  k  1  —  :  .  l  =  1 1  ]  • 

.ft. 

We  have  now 

(-hk2/hki)XkiUk  =  Xkl(-hk2/hki)l  T  =  {~hk2/hki)Hk]  =  -Hk  2. 

Lsing  the  block  submatrices  we  now  write  this  equation  and  the  equation 

AT  2  V*  =  Hk  2 
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-Yjt  11  (-hki/hujl’ki  +  Xk  12[  —  h k2 / h kl)[  k2  -  -Hk2 
Xk  21^*1  +  Xk  22  —  Hk 2- 


Adding  these  two  equations  we  obtain 


( —  h-k2  /  h  k\  )£l 

-'■Aril  -'A:  1 2  -'*22!  '*1  —  (^Atf/^Arl  K'A:2 

1 

(-hk2/hk\)(] 

Lk^\  =  Vk]  -  (h  k2  /  h  ki  )Ijt2 
1 


^Ar-t  I  1  -  !  Tk  +  1 


T2k~3  i  t:t+  1 


the  scalar  product  of  the  row  vector  [  -  2  •••  ^2*4-3]  and  the  column  vector  Uk+i- 

Adjoining  the  row  vector 


to  the  matrix 


we  have 


■ector  1 

Cjt-2  • 

•r2k->-3 

k*- 

2  •  •  ■ 

* 2k  4  3  ] 

[-v,„ 

Afc  )2 

A"  A:  22  ] 

•Cl 

^2 

1 

* 2 

•^3 

.r* .  2 

•CA :  *  2 

Tk^  3  ■ 

J2A  ■  2 

lAr4  11  — 


The  Eqs  (49)  and  (50)  show  the  computations  involved  in  going  from  the  kth  stage  to  the 
(k  +  1)**  stage  for  the  matrix  A\  j .  If  + 1  1  =0  we  are  finished,  if  not.  we  have  to  obtain 
the  vector  I’*-  j  and  the  value  hk  +  \  2  corresponding  to  the  matrix  AVi  2- 

Delete  the  first  row  of  the  matrix  AVij  and  denote  the  matrix  arising  from  this 
deletion  by  A**.  Set 

H  =  -Y*r*_i  (51) 

and  note  that 

Xk2(-hk,u/hk7)Vk  =  -H-  (52) 

IK 


Note  also  that 


where 


Now  we  can  write 


X*  =  |X*2  -V 1 


•  *lk  .  3  • 


r*+I  = 


and  adding  Eqs  (51)  and  (52)  we  obtain 


X*  \Uk  +  n_  (v>n/fcM)ni 


i/  _  f'*+i  i  -  (^t  +  i  i/hkz)} k 

Vk  + 1  -  j 


^*+12  ~  I  3  •••  *2*4  4  ]''*■+ 1- 


Adjoining  the  row  vector  [  .r*  +  3  . . .  x2*44  ]  to  the  matrix  X*  we  obtain  the  matrix  Ar*  + 1 2- 

Once  the  progressive  algorithm  is  started  the  essential  computations  are  those  given 
by  Eqs  (49)  and  (50)  and  Eqs  (53)  and  (54).  On  examining  these  equations  one  sees  that 
the  progressive  algorithm  could  fail  if  for  ft*  j  ^  0  the  value  ft*2  =0.  For  our  application  it 
is  clear,  in  principle,  at  least,  that  such  a  condition  cannot  arise,  since  for  an  appropriate 
value  of  k  both  X*  1  and  X*2  should  determine  the  coefficients  of  the  polynomial  P(t/).  Eq 


Now  it  is  clear  also  that  in  numerical  work  wre  cannot  expect  to  obtain  ft*  ]  =  0.  Thus 
in  practice  the  algorithm  will  terminate  when  the  value  ft*  j  is  sufficiently  small. 

We  will  now  summarize  the  above  algorithm,  using  superscripts  to  distinguish  one 
stage  of  the  computations  from  the  following  stage.  The  algorithm  is  as  follows: 
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For  A-  >  1 


tf”  =  -*a/*2  and  «,3f!l) +  »4  =  *j ’• 

di+i "  =  (-d*’/»!‘')d*'  +  <!*’ 


<1*”’  = 

dV,"  =  (-4,|/M,’i  +  <!*’ 

+ - t-  -r2i.2£[*Vt  + 


(  £  1  \ 

If  hj  J  is  numerically  sufficiently  small,  finished.  If  not,  then 

?r,)=er,)-^r,)Mt) 


(fc  +  l)_*(fc+l)  ,(*:) .  (*-l) /,  (*) 

—  ^k  H  nl  /  "2 

+ 1)  .(^1)/fc(*) 

>*+1  ~  ^*-*-1  "l  /"2 

and 

*4  ^  =  rk  +  3(]k  ^  + - V  r2k^3^\ ^  +  T2k^4-  (60) 

We  note  that  the  progressive  algorithm  is  equally  applicable  if  the  in  Eq  (37)  are 
replaced  by  the  averages  x}.  Eq  (34). 
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SECTION  IV 


THE  METHOD  OF  LEAST  SQUARES 

In  this  section  we  examine  in  some  detail  the  method  of  least  squares.  For  points 
x*  =  (xj *,...,*«*)  in  n-space  and  y*  =  {ylk . ymk)  in  m-space,  we  suppose 

y*  =  </(**)• 

The  function  #(x)  is  either  not  known  or  not  readily  evaluated.  For  some  simpler  function 

y  =  /(*•«) 

depending  on  a  parameter  set  a,  the  parameter  set  a  is  determined  so  that  the  sum  of  the 
squares  of  the  magnitude  of  the  error,  y *  -  y*  =  y*  -  f(xk,a)  is  least.  Because  of  our 
special  interest  this  review  only  examines  the  case  y  depending  linearly  on  x. 

For  this  case  we  regard  the  points  xk  and  y*  as  n-vectors  and  m-vectors  respectively. 


x*  =  |  xu 


*nt  r ' 


y*  =  (yu  Vmkf 


and.  for  example. 


lly*l!2  =  Y  y)k- 

We  will  limit  our  discussion  to  the  case  where  the  number  N  of  such  vectors  is  finite  but 
much  greater  than  either  m  or  n.  For  an  (m  x  ??)  matrix  A  set 


and  the  error 


y*  =  Axk 


=  y*  -  y*- 


Desired  is  a  matrix  A  for  which  the  average 

^[||e*i!2)-r[||y*-ytj|2] 
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is  minimal.  For  k  =  1,2, ...  we  obtain  from  Eqs  (61)  and  (62)  the  system  of  equations  for 
the  first  component  of  the  error 


elk  —  V\k  ~  °11  x\k  ~  a)2x2 k  ~  ~  alnxnk  (64) 

and  then 

-V  A  n  j 

E  4  =  E  ( y i*  -  E  a\jxjk)  ■  (65) 

k~i  t  =  l  j=l 

In  this  equation  we  have  a  function  of  the  parameters  au , . . . ,  a,„.  We  know  from  the 
calculus  that  to  find  an  extremum,  we  differentiate  with  respect  to  a j / .  equate  to  zero 

for  /  =  1 . n.  and  solve  the  resulting  system  of  equations  for  the  a j,.  We  have  on 

differentiating  Eq  (65)  with  respect  to  au,  for  1  <  /  <  n 

A  n 

Y^(yn*lk  -  'S2<*i}-r}kxlk)  =  °-  (66) 

fc=l  )=] 

These  equations  can  be  rewritten  as 

1  A  j  ^  Y  A 

(  ft’  E  xlk')al]  +  •  •  ‘  E  Xnk  Xlk)a\n  ~  ft  E  Wit  xlk-  (67) 


Bzz(jJ)  —  »T  E  Tjk  xtk  (68) 

‘  k  1 
1  \ 

Byz(lJ)  =  ft  E  S/itt  (69) 

k  -  \ 

The  functions  Bzz(j,l )  and  Byz(j.l )  are  called  the  correlation  function  and  the  cross¬ 
correlation  function  respectively.  When  there  is  no  danger  of  confusion  the  subscript  xx 
will  be  omitted. 

The  Eqs  (67)  determine  the  first  row  of  the  matrix  A  so  that  the  average  of  the  square 
of  the  first  component  of  the  error  is  minimal.  In  the  same  way  the  remaining  rows  of  A 
can  be  determined  to  minimize  the  average  of  the  square  of  the  corresponding  component 
of  the  error.  Using  the  notation  just  introduced  for  the  sums,  the  equation  for  the  desired 
matrix  A  can  be  written  as 

f0ll  ...  aIn  1  [  Bzz(l.l)  ...  BZ1  ( 1 ,  n )  1  fB„(l.l)  ...  B„I(l,ri)l 


a\  1 

®ln 

aml 

amn 

5„(n,l) 


B«(l.") 

Bzz(n.n ). 


■  Byz{m,  1)  ...  ByT(rn.n). 


no 


It  is  clear  from  Eq  (68)  that  the  coefficient  matrix  of  A  is  symmetric. 


It  is  convenient  at  this  point  to  obtain  a  relation  which  will  be  useful  later.  Let  us 
write 


Then 


elk  ~  (j/lfr  ^ ",  alixik)  (Vile  all^lt  a\nrnk) 

J  = ' 


=  y\k  -  XI  auxikV\k 


n 

-  All  (l/U* It  -  XI  aljxjk*\k) 


J=1 


II 

—  •  ••  —  fij  n  lik^nk  ~  XI  a\ ]x  )kxnk}  • 


t=l 


XI  e?t  =  XI  (j/u  -  XZ  fli>*>fcVi*) 

t=i  *=i  j  =  i 

A’  n 

-  au  XI  (yiktik  -  XI 


*- 1 

A  n 

- aln  XI  (vik*nk  -  XI  a\jrjk*»k)- 

*  1  j=  1 

It  follows  from  Eqs  (66)  and  (71)  that 

n 

E\e\k\  =  ^[yuyu]  -  XI 


(71) 


a 
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wav  the  vectors  y *  are  regarded  as  the  columns  of  an  (m  x  A’)  matrix  and  the  transpose 
of  this  matrix  is  V.  Let  X}  and  Y}  denote  the  jth  columns  of  A'  and  Y . respectively. 

Aj  =  [  Jji  Ij 2  ■  ■  ■  T}\  ] 

y}  =  [yj\  y} 2  •  ■•  y;A-]r- 

That  is  the  kth  component  of  X}  is  the  jth  component  of  x and  similarly,  the  kih  compo¬ 
nent  of  Yj  is  the  jth  component  of  y*.  The  vectors  X}  and  Y}  are  N-dimensional  vectors. 
A  scalar  product,  { Xj.Xi )  or  (Y},Xi).  for  these  vectors  is  given  by  Eqs  (68)and  (69).  It 
is  obvious  that  the  conditions  for  a  scalar  product  are  satisfied,  since  the  scalar  product 
defined  here  differs  from  the  customary  scalar  product  only  by  the  factor  1/AT.  Set 

l|A>|&  =  (A,,  A,).  (74) 

Consider  the  system  of  equations 


^  j  —  ajiXi  +  •  •  •  +  aJnX„.  (75) 

This  is  an  over  determined  system.  It  has  a  strict  solution  if  and  only  if  Y}  lies  in  the  space 

G  spanned  by  the  vectors  A'i - ,  A'„.  If  Y}  is  not  in  the  space  9  then  for  the  vector  Y}  in 

9  closest  to  Yr  the  difference  Y}  —  Y}  is  orthogonal  to  every  vector  A'/,  for  1  <  l  <  n. 

The  proof  of  this  assertion  is  by  contradiction.  Suppose  A'  is  a  vector  of  the  spanning 
set  A'), . . . ,  A'„  which  is  not  orthogonal  to  Y}  -  Y}.  Set 

b  =  (Yj-Y,.X) 

and 


(76) 


I  o  — ■  5  ,  +  •  —  —  A' . 
’  llA'Hm 


(77) 


Then 


\y}  -  yo\\2m  =  (y,  -  y,  - 


AM  '  -  I  '  - 


(h  ~YyY}-Yj) 

-  2b7  (yj  ~  Tj.A')  +  ----- 

!  VI 1 2  '  3  31 

% 

b2 


b  V  A') 
.’■112  ' 


A’llm 


(78) 


l!  A(|m 

\Y,  ~  Y,  II" 


All 2 


ilA'Hm 
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From  this  we  see  that  if  some  member  X  of  the  spanning  set  Aj .  ■  •  A'„  is  not  orthogonal 

to  Yj  —  Yj,  then  there  is  a  vector  V0  in  5  closer  to  Yj  than  Y}. 

Hence,  instead  of  trying  to  solve  the  over  determined  system.  Eq(75),  we  want  to 
determine  the  vector 

} j  =  ajiXj  +  -••  +  a}nX„  (79) 

closest  to  Y}.  We  have  then 

Y,  -  Yj  =  Yj  -  a)X  A, - ajnXn.  (80) 

Note,  for  j  =  1  this  equation  is  just  the  system,  Eq  (64),  above.  Since,  as  has  just  been 
shown,  this  difference  must  be  orthogonal  to  the  spanning  vectors  A'/  for  1  <  /  <  n,  we 
obtain  the  system  of  equations 

aji(XhXl)Y---  +  aJn(  X,,X„)  =  (XhY,).  (81) 

It  is  clear  that  we  can  rewrite  Eq  (81)  as 

’  BXI(IJ) 

[aj\  •••  G?n  ]  :  —  Byj(j.l). 

.  Bx,(nJ) . 

Doing  so  for  /  =  1 . n  and  j  —  1 . m  we  obtain  again  the  Eq  (70). 

A  difficulty,  albeit  minor,  with  this  last  derivation  is  that  it  focuses  attention  on  the 
vectors  X}  and  Y}.  There  is  a  tendency  to  lose  track  of  the  relation  of  the  result  obtained 

to  the  vectors  x*  and  y The  Eqs  (81)  determine  coefficients  ajX . ajn  and  in  turn  by 

Eq  (79)  a  Y}  for  which  ||Vj  -  Y}\\m  is  minimal.  The  kth  component  of  Y}  is  y]k,  hence 
denote  the  kth  component  of  Y}  by  y}k.  Then 

II  v,  -  y,|&  =  lSrY.(y^-yjk)2 

*  k  1 

m  .  X  m 

.7-1  •  fr-lj-l 

=  v  X  lly*  -  yt-ll2 

*  k=-  ) 

=  -Cflly*-  -  y*il2!- 


and 


From  Eq  (79)  we  have  for  the  kth  component  of  Y} 


l 

I 

j 

boo tJCUOUG 


Vjk  ~  a]lIlk  +  '  *  '  +  0]nxnk- 

Thus  for  j  =  1. . . . ,  m  the  system,  Eqs  (81)  determine  the  elements  of  a  matrix  A  for  which 


y*  =  Ax* 


and  E'flly*  —  y*||2]  is  minimal  as  desired. 
Again  from  Eqs  (61)  and  (62)  we  have 


e*  =  y*  -  Ax* 


=  \  Im  -A 


Post  multiplying  this  equation  with  [y*  x*  ]  we  obtain  a  matrix  equation  in  block  form 


fe,vr  eLxrI-f/  4lfy*y*  y*X* 
|e*y*  etx*l-l-'ru  -A  j  T  r 

xkJk  xkxk 


This  gives  a  system  of  matrix  equations 


«*y*  =  y*y*  -  Ax* y[ 
ekxl  =  y k*[  -  Ax*x* . 


For  the  left  side  of  Eq  (84)  we  have 


T 

e*x*  = 


ei  k*\k 


emkI\k 


ClkTnk 


Cmk*nk . 


(V,  -V„.Y,) 


Oi  -  V,..Vn) 


£[<****]  =  : 

.  0  m  “  l)  •  •  •  0  m  ~  ^ m-  An)  . 

Now  we  have  just  seen  that  for  the  Yj  closest  to  Y}  the  difference,  that  is.  the  least  error, 
must  be  orthogonal  to  the  spanning  vectors  A'j . A*„.  It  follows  that 

E[e*x[]  =  0  (85) 
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and  hence 


AE\xkx\\  =  E[ykxl\. 

That  is 


E\x\kxnk\ 

'  E\yxkxxk\  .. 

.  E\inkxxk\  • 

■■  E \xnkxtlk\. 

■  E\ymkxxk\  . 

■■  E\ymkxnk). 

and  this  is  just  Eq  (70)  again. 

From  Eqs  (87)  and  (70)  let  us  set 


E[xlkxlk]  .. 

•  E\xlkxnk\' 

BXI( 1,1)  .. 

■  Bxx(n.l)' 

E[xnkxllc] 

•  E[xnkxnk). 

Bxx(l.n)  .. 

■  Bxx(n.n). 

E\yikxxk} 

•  •  E{ylkxnk\' 

Byx  (  I '  1  ) 

Byx(l.n) 

E[ymkxlk]  . 

•  ■  E\ymkxnk] . 

.Byx(m.l)  . 

..  Byx(m.n) 

Then  we  have 

A  —  By  j  Bit* 

and  for  this  matrix  A  the  average 


(86) 


(87) 


(88) 


(89) 


T 

ykYk 

T 

y  k*k 

An  1 

T 

Lx*y* 

T 

x*x* 

l-^rJ 

BlIM2]  =  £(l|y*  -  y*ll2]  =  £!l(y*  -  Axk\\2 

is  minimal. 

From  Eq  (82)we  obtain 

—  [  An  —  A 

=  ly>y[  -  A*kvl  y -  Axkx\ ,  ^  _ 4, 

It  follows  from  Eqs  (86)  and  (83)  that 

E\ekek  \  =  E\ykyl  -  AE[xkyl\ 

=■  Byy  -  A  Bj  y  =  £[e*y[]. 

Now  from  Eq  (73)  the  ith  diagonal  element  of  the  right  hand  side  of  this  equation  is  the 
average  of  the  square  of  the  ith  component  of  the  error.  Hence  we  have 


(90) 


trace  (By  y  -  .4B7y)  =  || 2]  - 


(91) 


There  is  a  difficulty  with  the  method  of  least  squares  which  we  will  now  try  to  show. 
For  any  n- vector  a 

(arx*)2  =  a7x/txta  >  0.  (92 J 

From  this  it  is  clear  that  the  matrix 

B/t  =  £>***)  (93) 

is  symmetric  and  nonnegative.  Hence  there  are  eigenvalues  A  i  >  A2  >  . . .  >  A„  >  0  and 
associated  eigenvectors  Uj,. . .  ,u„  satisfying  the  conditions 

B„u2  =  A;iij 
u,rUj  =  St} 

(94) 

=  1  for  i  —  j 
—  0  otherwise. 

Then  we  can  write 

n 

Bji  =  Y.  A*u*u*  (95) 

k  1 

and  if  A„  >  0  then 

^  1 

=  E  AtUfcU*'  (96) 

Consider  the  product  BTVa.  Since  in  practice  the  vector  a  is  experimentally  determined 
there  is  error  in  a.  If  now,  for  example,  An  is  much  less  than  1  then  the  error  in  the 
component  of  a  in  the  direction  u„  will  be  greatly  magnified.  A  similar  remark  holds  for 
all  values  of  k  for  which  A*  is  much  less  that  I. 


’  it'  , 
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SECTION  V 


LATTICE  FILTERS 


In  the  previous  section  we  discussed  the  method  of  least  squares.  In  this  section  we 
will  see  the  simplifications  to  this  method  when  the  data  is  stationary.  Lot  x(t)  denote 
a  real  stationary  random  sequence  then  the  mean  is  m(t)  =  £'|x(/)]  and  for  a  stationary 


process 


m{t  +  r)  =  m(t)  for  all  r. 


From  this  it  follows  that  m(t)  is  a  constant  and  it  is  sufficient  to  determine  m(t)  for  t  -  0. 
In  the  following  we  shall  suppose  for  the  sequence  x(t)  that  the  mean  is  zero  for  if  not,  we 
can  consider  instead  the  stationary  random  sequence 

*«)-£[*(')!• 

Thus,  here  and  in  the  following  we  have  an  ensemble  or  set  of  sequences  and  x(t)  is  a 
representative  member  of  this  set.  The  mean,  rn(t)  is  the  average  over  the  sequences  in 
this  set  of  sequences  at  a  fixed  value  of  t. 

The  correlation  function  D(f.s)  i« 

B(t.s)  =  E[x(/)x(s)]  =  E[x(s)x(/)]  =  B(s,f)  (98) 

and  for  a  stationary  process 

B(f  +  t.s  4-  r)  =  B(t.s)  for  all  r.  (99) 

Taking,  for  example,  r  =  -s  we  see  that  B(t.s)  depends  on  the  difference  t  -  a  alone. 
Hence  for  a  stationary  process  we  can  write 

B(t,s)  -  B(t  -  a,0)  =  B(t)  where  r  =  t  -  s.  (100) 


The  correlation  function  has  the  properties: 


5(0)  >  0 
B(-t)  =  B(t) 
|5(r)|  <  5(0) 


(lOlo) 

(1016) 

(101c) 
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Properties  (a)  and  (b)  follow  immediately  from  Eqs  (98-100)  and  (c)  follows  from 

E[{x(t  +  r)  ±  -r(/))2]  =  E\xl(t  +  r)J  1  2 E\s(t  +  r)x(t)}  +  E\x2(t)}  >  0. 

Lot  us  note  also  that  for  any  n.  a, - ,  a„.  1 1 . ln  that 

(o,j-(/,)  +  ••  -  +  a„x(t„))2 

*(Mi  for 

=  [«i  On]  :  |*(*i)  ••  *(<n)j  • 

.x(<„)J  U„. 

This  equation  can  be  written  as 

(o13'(/1)  +  •  •  •  +  a„x(t„))2 

rx(/.1)x(r,)  ...  x(<, )x(<„) i  ra,' 

=  [fli  •••  On]  i  :  : 

.jr(/„)x(/,)  *(<nM*n)J  Un. 

Taking  the  expected  value  of  both  sides  we  obtain 

•  B(0)  ...  B{t,  —  <„)  I  o, 

(«,  ...  a„]  :  •.  i  i  >0.  (102) 

.B(tn  -  /,)  ...  5(0) 

We  shall  refer  to  the  matrix  of  this  quadratic  form  as  the  correlation  matrix.  From  the 
properties  of  the  correlation  function  B(t)  it  is  clear  that  the  correlation  matrix  for  a 

stationary  process  is  symmetric.  Usually  t, . t„  are  successive  values,  for  example, 

/  +  1 . t  +  n.  and  the  correlation  matrix  becomes 

’  B{  0)  ...  B(1  -  n ) ' 

(103) 

.  5(n  —  l)  ...  5(0)  . 

Note  that  the  element  in  the  j'h  row  and  kth  column  is 

B{t  +  j.1  k)  =  B{j  -  k) 

and  the  element  in  the  (j  +  1)J/|  row  and  the  ( k  +  l)(/l  column  is 

B(t  +  j  +  1./  +  k  +  1)  =  B(j  -  k). 
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From  these  observations  it  follows  that  the  correlat  ion  matrix  determined  from  a  stationary 
random  sequence  is  a  symmetric,  nonnegative  Toeplitz  matrix.  In  this  case  the  correlation 
matrix  is  known  if  one  knows  the  elements  in  the  first  row  or  column. 


Consider  a  scalar  stationary  random  sequence  x(t)  with  zero  mean.  Set 

T{t±})  =  X(t±j). 

Thus  . .  „_!) . Xp_  j),ijt),X(t+1), . . .  is  a  set  of  scalar  values.  A  linear 

forward  predictor  of  order  n,  LFP(n),  is  defined  by  the  equation 

xt  =  anix(t-l)  +  an2x(t-2)  +  •••  +  annx[t-n)-  (104) 

Here  xt  denotes  the  predicted  value  of  x<  based  on  the  n  values  xf_j, . . .  ,x<_„.  Let  eni 
stand  for  the  prediction  error  of  order  n. 

n 

e„(  —  X(  —  Xf  =  X(  —  anjXt~y  (105) 

j=i 

The  prediction  coefficients  «nl .....  a„„  are  determined  so  as  to  minimize  the  prediction 
error  in  the  least  squares  sense. 

To  relate  this  to  our  previous  discussion  of  the  method  of  least  squares  set 

X  —  [  Xj  _  I  ...  X(_n  j  . 


This  n-dimensional  vector  x,  which  is  a  representative  of  many  such  n-vectors,  plays  the 
role  of  the  n-vectors  x*.  The  representative  value  x*  plays  the  role  of  the  m- vectors  y*  and 
the  n-vector 


the  role  of  the  matrix  A.  Having  noted  this  connection  to  our  previous  discussions  we 
write  Eq  (105)  as 


ent  ~  xt  xt  —  [  1  an I 


=  11 


x1-n  J 


(106) 

(107) 


3] 


Post  multiply  Eq  (107)  with  the  row  vector  [  x,  ...  xt  „  ]  =  [  xt  xT  )  and  obtain,  as  in 
the  previous  section. 


[<-ntXi  en<xrj  =  (l  -a], 


XtT,  x,xr 
X  X,  XXT 


(108) 


Taking  the  expectation  we  know  from  our  discussion  in  SECTION  IV  that  we  can  obtain 
the  two  equations 


Ejxfx,]  -  a lE\xxt\  =  E\tnixt\  (109) 

E\xtxT}  -  al  E[xxt]  =  E[entxT\  =  0.  (110) 


Eq  (110)  determines  the  vector  a„  and  once  this  vector  is  known  Eq  (109)  determines  the 
value  E\entxt\. 

From  Eq  (107)  we  obtain 


^nt^nt  —  [  1 


x,x,  x,xTl 

1 

.  XX(  XXT  \ 

~«n 

(in) 


Taking  the  expectation  and  performing  the  indicated  multiplication  on  the  left  we  have 
E[e2nt\  -  \  E[x'{]  -  aj,  r[xx,]  E\x,xT}  -  a lE[xxT] 

It  follows  from  Eqs  (110)  and  (109)  that  this  equation  becomes 


1 


E\e2„t}  =  E[x(2]  -  alE\xxt\  =  E\enixt\. 


(112) 


Set 


rl  =  EK). 


(113) 


Finally,  then,  taking  the  expectation  of  Eq  (108)  we  can  write  the  two  equations,  Eqs  (109) 
and  (110)  as  a  single  matrix  equation 


\r2n  0  ...  0  j  =  [  1 


5(0) 

B(-l) 


5(1) 

5(0) 


B(  —  n)  5(1  -  7 ?) 


5(h) 
5(h  -  1) 

m  . 


(114) 
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Let  us  denote  the  correlation  matrix  in  Eq  (114)  by  B„^  j.  We  observe  that  except  for 


a  normalization  factor  the  row  vector  !  1  —a, 


-u„„  )  is  the  first  row  of  the  inverse 


of  Bf,-).  Also  from  examination  of  Eq  (114)  we  see  that 


rl  0  ...  01 

0  ...  0  rl  J 


Our  object  is  to  show  that  if  we  have  the  prediction  coefficients  ani . antl  of  order  n 

then  the  prediction  coefficients  of  order  n  -4-  1  are  readily  obtained. 


We  suppose  the  correlation  matrix  B„  is  known  for  every  value  n  of  interest.  We 

suppose  the  prediction  coefficients  an] _ ,an„  are  also  known.  Then  from  Eq  (115)  the 

value  rl  is  known  and  we  can  compute 


1  ®nl 
0 


0fin  ^  11  _ 

1  Un +2  ~ 

an  1 


_  f  rn  0  ...  0  qn 


qn  0  ...  0  rl 


The  right  side  of  Eq  (H6)  is  not  quite  of  the  right  form.  There  is  a  value  q„  which  would 
be  zero  if  we  had  the  right  prediction  coefficients  of  order  n  +  1.  This  value  of  q„  is  given 


<7n  —  D(u  +  1)  'jr  dnlcB(rt  i  1  k) 
i 


It  is  clear,  however,  that  if  we  multiply  Eq  (116)  by  the  matrix 


-9n  ’rl 


Qn/rr)  1 


the  right  hand  side  takes  on  the  right  form.  The  prediction  coefficients  of  order  n  +  1  are 


given  by 


®n+  1  j  —  ®  n  j  (Qn  /  rn)an  n-r  I  -  j  ft>r  1  5;  j <  H 


*  I  n^l  —  Qn/rn- 


One  has  also 


.2  _  _2  _2  /  2 

rn  -  1  ~  rn  <?»i  /  rn- 


(120) 


Starting  values  are  readily  obtained.  We  have 


1  — O] ,  j  [  B( 0)  £(1)1  _  rf  0 

-ou  1  l[B(-l)  0(0),  0  r] 


s 


B? 


& 

1^.1 


>*i  i«c^f:i 


»**  -ti,.*i,-t»,.1fl  ■>.» li’ .  >»»  _*. 


ft 

a 


f 


From  this  equation  we  obtain 


-«,,£(())  +  Z?(I)  =  0 


ri  =  -0(0)  ~  «ii-0(l) 

_  D'l( 0)  -  B2(l)  , 


Consider  now  a  stationary  random  sequence  x(/)  of  m-vectors.  A  linear  forward  pre¬ 
dictor  of  order  p  is  of  the  form 

x(t)  =  Aplx(t  -  1)  +  • •  •  +  Appx(t  -  p) 

fx(':1|l  (i24) 

—  ( Ap\  . . .  App  ] 

.x(l  -  p). 

Here  the  lengthy  mp-vector  [x7^/  -  1)  ...  xT(t  -  p)  ]7  plays  the  role  of  the  n-vector  x* 

of  Section  IV.  The  (m  x  rnp)  matrix  [  .4pl  . . .  App]  plays  the  role  of  the  matrix  A  and 
x(M  the  role  of  the  vector  y*.  Again,  for  the  stationary  case  the  mean 

£'[x(l)]  =  E\x(l  +  r)]  for  all  r. 

Hence  as  for  the  scalar  case  we  may  suppose  x(/)  has  mean  zero. 


The  error 


«*P(0  =  x(0  -  *(0 


““  [An  -^f)l  Appl  ’ 

U(i  -  p)  J 

Post  multiplying  this  equation  by  the  vector  [  x^(/)  ...  xT {t  -  p)  ]  obtain 

[ep(0xr(/)  ep(t  )xT  (t  —  1 )  ...  ep(t  )xT  (t  -  p)  ] 

for  the  left  side  and  the  matrix 

x(/)xr(t)  x(/)xr(/ -  1)  ...  x(t)xT(t-p) 

x(l  -  l)xr(0  x(l  -  i)xT(p  -  1 )  ...  x(l  -  1  )xT(1  -  p) 

x(/-p)xT(/)  xO-p)xT(t-  1)  ...  x(t  -  p)xT (t  -  p) 


L__ 


multiplies  the  matrix  [  /,„  -Ap\ 


-App]  on  the  right.  Set 


Iip(l)  =  £[x(/)x7V)]. 


(125) 


Taking  the  expectation  of  the  equation  just  described  we  obtain 

[tfp(0  0  01  = 


I  m  A  pi 


A 


ppi 


Bm(  0)  Sm(l)  ...  Bm(p) 

Bm(-l)  Bm( 0)  ...  Bm(p  -  1) 


(126) 


. Bm(-p)  Bm(l  -  p)  ...  £m(0) 

The  subscript,  m,  is  to  indicate  that  the  entries  Bm(j)  are  matrices  .  This  matrix  in  this 
block  form  is  a  block  symmetric  Toeplitz  matrix  and  a  procedure  like  that  described  by 
Eqs  (115-118)  can  be  used  to  determine  the  prediction  coefficients  Ap\,. . .  ,App 
Starting  values  are  readily  obtained.  We  have 


A})  =  Bm(l)Bm  (0) 


(127) 


and 


*,(0  =  Bm(0)-,4nflm(-l).  (128) 

Set  Bmp+i  equal  to  the  matrix  with  block  elements  Bm(j).  Eq  (126).  Then  as  in  Eq  (116) 

(129) 


*m  Ap  | 

-App  0 

B  ,  - 

RP(t)  o  .. 

o  Qp 

.  0  -App  .. 

—  Ap l  It n  . 

“rn  p->  2  — 

QP  o  . 

•  0  /?„(/)  J 

From  this  equation  we  have 


QP  =  Bm(p  +  1)  -  Y,  ApkBm(p  +  1  -  k). 

k  I 


(130) 


Multiply  Eq  (129)  on  the  left  by 


K, 


-QpBp'iO 


-QpRp'V) 


(131) 


1 


& 


■  D.  l*i  •’i.t'iJt.l 


1 

ivJ 

I 

$ 


for  I  <  j  <  p  and 


-■f  ;j  <  i  p  *  i  —  QPRp 1  (0- 


In  practice,  since  Rp{t)  and  Qp  are  known  one  first  solves 


Api  i  p *  i  Rp(t )  —  Qji 


for  4p^]  p_) .  One  then  uses  this  matrix  in  Eq  (132)  to  compute  the  matrices  .4p+)j 

Let  us  review  first  the  scalar  case.  The  attractive  features  of  the  stationary  process 
are  that  one  does  not  have  to  choose  in  advance  the  order  of  the  system  for  determining 
the  predictor  coefficients.  It  is  not  necessary  to  compute,  store  and  solve  a  system  of 

equations  involving  the  matrix  B„.  It  is  sufficient  to  know  the  values  B( 0) . B(n)  of 

the  correlation  function  for  a  sufficiently  large  value  of  n.  As  the  prediction  coefficients 
a„\, ... .  ann  are  computed  for  n  —  1.2....  one  is  able  to  monitor  the  results  by  means  of 
the  value  r£.  The  value  r 1  is  the  average  of  the  square  of  the  error.  This  value  can  be 
used  as  a  criteria  for  determining  the  order  of  the  linear  forward  predictor.  If  for  some 
value  of  n  one  obtains  <  0  we  infer  that  the  experimental  data  does  not  support  the 
determination  of  prediction  coefficients  for  a  value  of  n  this  large. 

Corresponding  remarks  hold  for  the  vector  case.  Here  again  we  do  not  need  the  matrix 

B„,p  but  its  entries  Bm(0) . Bm(p)  for  p  sufficiently  large.  The  process  of  determining 

the  prediction  coefficient  matrices  .4p] . App  for  p  =  1.2....  is  monitored  by  means  of 

the  diagonal  entries  of  the  matrix  Rp(t). 


.vyjvwv 
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SECTION  VI 

CANONICAL  VARIATE  ANALYSIS 

Here  we  review  a  method,  canonical  variate  analysis ,  for  determining  a  one  step  ahead 
linear  forward  predictor  which  is  optimal  with  respect  to  a  prescribed  norm.  As  above, 

SECTION  IV.  we  are  concerned  with  a  set  of  n-vectors  Xj, _ x jv  and  a  set  of  m-vectors, 

yi . y,\  ■  The  method  does  not  require  that  the  x*  and  y*  be  stationary  and  jointly 

stationary  as  in  the  previous  section.  However,  stationarity,  if  present,  will  provide  some 
simplifications. 

Let  G  denote  a  symmetric  positive  definite  matrix  of  order  m.  Then 

llytllc  =  y*G’y* 

=  y  lG'l*G'i*yk  (135) 

=  l|G,/2y*||2. 

Observe  for  any  orthogonal  matrix  V  of  order  m  and  for 

a*  =  V tG'  2yk  (136) 

that 

!lz*!!2  =  l|y*llc-  (137) 

Let  us  suppose  that  the  matrix 

=  £[x*x*]  (138) 


is  positive  definite.  Then  is  nonsingular.  Here  we  want  to  determine  a  matrix  J  of 
order  n  which  transforms  the  vectors  x*.  into  vectors  w*. 


w*  =  ./x*. 


(139) 


for  which  the  covariance 


C(^.wt)  =  Bu.u,  =  /„. 
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(140) 


W  TTUirUTTUH 


^’(w*-wt)  =  £’lwi'w[] 

=  E[JxkxJjT]  (141) 

=  jbJ7jt. 

Since  the  matrix  BJT  is  symmetric  and  positive  definite  our  objective.  Eq  (140).  will  be 
realized,  if  for  any  orthogonal  matrix  V  we  take 

J  =  VT  B„1/2.  (142) 


The  covariance 

C{wk,zk)  =  E\wkzl]  =  BU)Z 

=  £(V'TB,V/*xfcyfG,/2r7)  (143) 

=  l,TV«l/2B,„G,/2r. 

If,  now,  U  and  V  are  the  orthogonal  matrices  in  the  singular  value  decomposition  of  the 
matrix  Bit*  2BJ!/(V1/2  then 

V'rB„,/2&7S,<V,/2r  =  S  (144) 

where  S  is  an  ( //  x  m)  matrix.  All  the  entries  in  5  are  zeros  except  for  the  diagonal 
elements  for  1  <  j  <  r.  where  r  is  the  rank  of  the  matrix  BJV.  and  r  <  l  ~  mtn[m.n). 
These  diagonal  elements  satisfy  the  inequalities  an  >  ,s22  >  •  •  •  >  srr  >  0. 

We  now  want  the  matrix  A  for  which 

i-k  =  -4w  * 

. f  l.e  1  (145) 


ck  =  7.k  -  ik  I  Im  -A] 

and  E  [  ||c*.  j |2  j  is  least.  As  we  know,  from  what  has  been  said  above,  taking  the  average  of 

[<****  1  =  [  Im  -A]  ®***r  Z*W*T 

’wk7.k 

obtain  the  system  of  equations 

B--  -  AB„  ;  =  £  [<***!] 

-  ABU,„,  -  E[etw[| 
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For  the  desired  matrix  A,  Ele^wJ J  =  0  and 


.4BU,U.  =  B.u, 


(146) 


and  since  Bu.„  =  Iri 


A  =  B,„,  =  ST 


by  Eq  (144).  Thus  the  predicted  value  i k  is  given  by 

ik  =  ST  wt 

r 

=  Y.snwjk 


and  one  solves  the  equation 

for  the  predicted  value  yk. 

Also,  as  we  have  already  seen, 


UTG'!*yk=ik 


—  I  T 

A  1 

T 

Tkz[ 

i 

7- 

l  1m 

A\ 

T 

* 

*- 

$ 

-  A7"  j 

From  this  equation  we  have 

^(e*e*]  =  (B„-  ABU(Z  B  -  A  Btt,u,  J 


=  B,.  -  AB„ 


because  of  Eq  (146).  Now 


and  then 


AKz  =  Biu,B,„,  =  SST 

Trace  E\ekeTk  \  =  £[|je*||2] 

=  T race  B-2  —  T race  SST 


=  T race  B--  -  £  ,s;2; 

}  i 


(147) 


(148) 


(149) 


(150) 


(151) 


(152) 


We  want  to  examine  the  observation.  Eq  (13).  in  greater  detail  by  considering  the 
scalar  case.  Application  of  the  method  to  the  scalar  case  is  relatively  straight  forward. 
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Thus,  for  a  random  sequence  with  zero  mean  and  for  a  provisional  value  of  n  we  want  to 


determine  values  otl] . ann  for  which  the  predicted  value 


xl  —  an  1  •*"(<  -  ] )  "t"  (ln'2x(t  2)  d"  '  ’  ’  "h  an rj  -r ( t  -  n) 


is  best. 

For  a  value  n 


Brr  =  E 


x(/-I)x(f-l)  ••  x(t-l)x(i-n) 

{.*(»- w)T(t-I)  ••  x(t  n)x(t-n) 


and  we  can  write 

B„  =  + - 1-  Af,unii£ 

where  A]  >  A2  >  •  •  •  >  A„  >  0.  Here  we  want  a  rather  large  value  of  n  for  which  the 
smallest  eigenvalue  An  >  0  and  not  too  small.  For  such  a  value  of  n  the  matrix  B„  is 
nonsingular  of  order  n.  The  matrix 

x{t  i)xt 

=  E  : 

.  x(i  -n)xl  . 

in  this  case  is  an  (»  x  1)  matrix  or  an  n- vector.  For  the  scalar  case  the  matrix  G  is  a 
positive  constant  which  for  our  purposes  here  we  take  as  1. 

We  readily  see  that  the  singular  value  decomposition  of  the  matrix  Bit  Bry  is 

VThIJl2biyU  =  |.s„  0  ...  of 

where  U  =  1  and  hence  the  A  matrix  is  the  row  vector 


A  =  {su  0  ...  0], 


(153) 


Here  one  must  keep  in  mind  that  in  the  present  case  t,  plays  the  role  of  the  vectors  y *  =  * * 
and  x*  =  [^p-i)  •••  •r(/~n)f-  Doing  so.  it  is  clear  then  that 


From  this  equation  one  could  obtain  the  coefficients  anl , . . . ,  ann.  On  the  other  hand,  if 
we  set 

x(t  IJ 

wt  =  \'T^>ii  1  '■  (155) 

■  X«  n). 

then 

Zt=«ll»'\k  (156) 

that  is,  the  product  of  and  the  first  component  of  w*. 

In  review  and  conclusion  then  the  main  feature  of  this  report  is  the  application  of 
Lanczos'  Progressive  algorithm  in  the  method  proposed  for  the  determination  of  complex 
exponentials.  From  our  limited  numerical  experimentation  with  the  method  we  believe  the 
method  merits  further  investigation  and  testing. 

The  method  of  lattice  filters  for  prediction  in  stationary  sequences  is  attractive  because 
of  the  simplicity  of  the  computations  and  the  ability  to  monitor  how  good  the  predicted 
value  is  by  means  of  the  quantity  r£  or  Rp(t).  The  matrix  square  roots  and  the  singular 
value  decomposition  in  the  canonical  variate  analysis  method  indicates  that  the  computa¬ 
tions  in  this  method  are  more  extensive. 
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